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GENERATION OF 

THREE DIMENSIONAL BODY FITTED COORDINATES 
USING HYPERBOLIC PARTIAL DIFFERENTIAL EQUATIONS 

Joseph L Steger and Y. M. Rizk 
SUMMARY 

An efficient numerical mesh generation scheme capable of creating orthog- 
onal or nearly orthogonal grids about moderately complex three dimensional con- 
figurations is described The mesh is obtained by marching outward from a user 
specified grid on the the body surface. Using a spherical grid topology, grids have 
been generated about full span rectangular wings and a simplified space shuttle 
orbiter. 


INTRODUCTION 


Body conforming curvilinear grids are often used m finite difference flow 
field simulations One reason for this is that the application of boundary conditions 
can be simplified in finite difference calculations if grid lines coincide with boundary 
lines. This is especially important in high Reynolds number viscous flow simulation 
in which high flow gradients near the body surface must be resolved 

The task of generating a satisfactory body conforming coordinate system is 
not easy. The grids must not be too distorted, they should have smooth variation, 
and they should be clustered to flow field action regions - typically near boundary 
surfaces. Moreover, the grids should be generated in an automatic manner that 
requires a minimum of user input. 

One approach for generating body conforming grids with minimum user 
input has been to solve a set of partial differential equations. In this technique 
level lines of £(x,y,z), rj(x,y,z), and f(x, y, z) that have monotone variation are 
sought as a solution of a set of partial differential equations. Generally, values of 
77 , and f are user-specified on the boundary surface, and constraints expressed 
as differential equations are used to develop the grid away from the boundaries. 
The most popular such approach requires the solution of a set of elliptic equations 
that satisfy the maximum principle ,1-5]; however, hyperbolic [6,7] and parabolic [8] 
governing equations have been used as well, at least in two dimensional applications. 

In this paper one way of extending the hyperbolic grid generation method 
of Steger and Chaussee [6i to three dimensions is developed. In two dimensions this 
grid generation method requires the solution of two partial differential constraints 


£xVx U £y*ly ® 


(la) 



txJIv - ZyVx = (AT) 1 (lb) 

where the subscripts x and y denote partial derivatives These equations are trans- 
formed to £. 77 computational space as 

J-In - ViVr = 0 (2a) 

x«j/r> — = A I (2b) 

and are solved by marching m 77 from an initial data plane 77(1. y) = const The first 
equation is a constraint of orthogonality. The second equation controls the mesh 
spacing with the user specifying the mesh control volume AT (actually area in two 
dimensions). A linearized version of Eqs. (2) is readily shown to be hyperbolic and 
suitable for marching in 77. Equations (2) are solved in computational space to give 
the x,y location of the specified £ = const and 77 = const grid lines 

The two partial differential equations, expressed as either Eqs. (l) or (2). 
have been referred to as a cell-volume procedure for grid generation For two di- 
mensional external flow simulation m which the outer flow boundary need not be 
precisely located, this cell-volume hyperbolic partial differential grid generation pro- 
cedure has been found to be an efficient way of generating orthogonal or nearly or- 
thogonal grids In the next section, a three dimensional extension of this procedure 
is developed 


THREE DIMENSIONAL GRID GENERATION EQUATIONS 

A body fitted exterior grid about an arbitrary closed boundary surface is 
desired. Only a simple topology such as that illustrated in Fig. (l) will be considered 
here. The body surface is chosen to coincide with f(x,y, z) = 0 and the surface 
grid-line distributions of £ = const and 77 = const are user-specified. The outer 
boundary f(x,y. 2) = f mox is not specified; it is only required to be sufficiently 
far removed from the inner boundary. Using f as the marching direction, partial 
differential equations are sought that produce surfaces of constant f, 77, and <r to 
form a nonsingular mesh system 

An extension of the cell-volume procedure to three dimensions is proposed. 
In three dimensions, however, there are three orthogonality relations and one cell- 
volume constraint. At any point four equations are available to predict the three 
unknowns x, y. and 2, so one equation must be discarded. Because f is the marching 
direction, it is natural to use only the two orthogonality relations that involve f 
this leads to the governing equations 

x f x ? -t ycy< + z^Z( = 0 

Xr,X<; -r yr?y ? + Z^Z, = 0 
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( 3 a) 

( 3 b) 



xty v Zr -t x^y^Zr, + x v y^Zc - x^y^Zfj - x v ycz : - x^y^zc = AV (3c) 
or. with f defined as ( x,y,z) f 

Fe-r c = 0 (4a) 

(4b) 

J~' = AV (4c) 


r n • r, = 0 


'd(x.y.z) 


•d(i i;.f) 

The first two equations represent orthogonality relations between f and f and be- 
tween 77 and and the last equation is the volume or finite Jacobian constraint. 

Equations (3) comprise a system of nonlinear partial differential equations 
in which x, y, and z are specified as initial data at f = 0. As developed below, 
linearization and analysis of Eqs. (3) about a nearby known state reveals that the 
system is hyperbolic with c taken as the marching direction. 

Let x°, y°. z° represent a nearby known state so that 

x = x" — (x — X°) = X° — X 

y = y u - y 

z = z° + i 

where x, y, and z are small. Substitution of these expressions into Eqs. (3) and 
elimination of products of tilde terms results in the locally linearized system 


Aq{t- r 0 )c -r B 0 {r - r 0 )r, + C 0 (r - fo) f = / 


(5) 


with 


and 
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Equation (5) is rewritten as 


Aorc - Bor „ — Cor. = e (7) 

with e = (0.0. AV — 2AV 0 ) < Now C~ 1 exists unless (AVu) — » 0 r. w’hich we will 
not impose, so (7) can be rewritten as 

C ~ 1 ^ u r £ - Cr l B.<r r , - f = C, -1 e (8) 

Although the algebraic verification is nontrhial and C ~ 1 B>> are symmetric 

matrices and are given in the Appendix (The matrices and the verification of 
symmetry were carried out using MACSYMA by Dennis Jesperson of the Ames 
Research Center.) Because C^ 1 Ao and Cq 1 Bq are symmetric, the linearized system 
Eq. (8) is hyperbolic and can be marched with f serving as the ‘‘time-like” direction. 

It can be pointed out that an analysis was attempted for the three orthogo- 
nality relations alone. These equations, however, are found to be improperly posed 
for marching with initial data in c Indeed, as best as we can discern, these three re- 
lations do not lend themselves to unique solutions regardless of the type of boundary 
conditions specified 


SOLUTION PROCEDURE 

The nonlinear system of grid generation equations given by Eqs. (3) are 
solved with a noniterative implicit finite difference scheme. An unconditionally 
stable implicit scheme is chosen so the marching step size m f can be arbitrarily 
selected based only on considerations of accurately generating the grid. Iterative 
solution of the nonlinear grid generation equations is avoided by expanding the 
equations about the previous marching step. As a consequence Eq. (7) is solved 
with the nearby known state 0 taken from the previous f step. 

The following subsections describe the numerical method for marching in f 
and describe the procedures for chosing the mesh-cell increments. Special differenc- 
ing must also be incorporated at axis singularities and at £ and r) boundaries. This 
detail is also described. 

Numerical Method 

Let A ^ = Arj = Ac = 1 such that £ = j — 1, n = k — 1, and ( — l — 1. 
Central spatial differencing of Eqs. (5) in £ and r\ with two-point backward implicit 
differencing in c leads to 

Ai6t[fi ^. | - r t ) -r BiS v {r l+ \ - fj) 4- CiV^r i+1 = g l+1 (9) 


where 


9l+i = 


( 0 
0 

^ AVj-r! j 
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and 


6tf 
s J 





[k-i - r fc _ i 
o 


V r r/_] = r,_! - r; 


Note that Aider i and Bid^fi were subtracted from the left-hand side of (7) to put 
the equations m a ‘delta* form Throughout onlv those indices that change are 
indicated, thus r/_] => r,.^;_jand j => r ; *j k /. etc 
Multiplying through by Cf 1 gives 


1 j4/^(F{ + i - r/) + C, - r t ) - 7(f/+i - r/) = C t Vt-ri (10) 

where 7 is the identity matrix. To reduce the inversion cost the difference equations 
are approximately factored as 


(I - Cf 1 Bi6 ri )(I - - U) = Cf'g l+ , (11) 


so that r; T i is obtained by solving sequences of one-dimensional-like block tridiag- 
onal systems 


(7 - C t 1 BiSr^g’ i + 1 = C { Vf+i 

(12a) 

(7 - Cf 1 Ai6c)V.ri^. 1 = g‘i+i 

(12b) 

ri ~ i = r t - V -r/_. 

(12c) 


In practice, numerical dissipation terms are added in the f- and r\~ directions. 
Typically, we have used a combination of fourth and second differences which are 
explicitly and implicitly included in the basic algorithm as 


J + Cr‘B,S n - t.,(AV),]|/ + Cr’A,{< - t, 4 (AV) { ](r l+1 - r,) = 

Cf ’sh , - lt. 4 (AV)f + q,(AV);|F, 1 

where, for example, 

(AV)„r = r k+ j - 2 r k -+- r* fc _i 

and 

(AV)?r = r^ 2 ~ 4 j- 6F, - Ar^-i + F,_ 2 

The dissipation coefficients are scaled to the local mesh using t e ^ — 0.5|!C _1 A|i, e er) = 
0.5|iC -1 7?j' and f, is 3 x t e . For simplicity, the matrix norms |jC~M|! and ||C -1 Bji 
have been approximated by a pseudonorm consisting of the square root of the sum of 
the squared elements. As an alternative to fourth differences, numerical dissipation 
terms such as 

e|AVr/|A Vr) 
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have been tried, but these are less effective than the fourth difference terms 

Following Kinsey and Barth 19 . additional smoothing and implicitness are 
put into the algorithm by differencing V-r = F as r/_j - r/ = (1 ~6)Fi~ \ -9Fi. with 
6 > 0 Ignoring the adjustments this entails in the user specified volume increments, 
which need not be exact, this differencing is e is incorporated Eq (13) as 


1 ~ (1 - 6)Cf l Bi6„ - AV)„ — (1 - 9)Cr'A,f>z - fl . c (AV)*!(r%, - r , ) = 

CT’ff-i - i < - e «( AV )? - £ e„(AV)‘in 

(14) 

Values of 0 of 1 to 5 are effective in preventing numerical breakdown in the case of 
highly concave body shapes. If some loss of grid orthogonality can be tolerated, 0 
can be set to 3 for most applications 

The coefficient matrices A\. B\ and C; contain £- and ^-derivatives which 
are formed using central differences. These matrices also contain derivatives for x c . 
y ( . and z- which are obtained from Eqs (3) m terms of £- and rj-derivatives. That 
is. Eqs. (3) are linear in the unknowns x-. y : . and c- The\ are easily solved as 


( X- ^ 

AV 

1 y^r, - y v zt } 

y- 

~ ( DetC ) 

X v z c — XcZri 

\ Z ‘J 



with 

Det(C) = [ycz n - y v zc) 2 + (x n z ( - xcz n ) 2 -+- (**y„ - x„y<) 2 

Note that (AV) 2 /(x? + y 2 + z 2 ) — Det(C) so that Det(C) will be zero if and only 
if the user specified AV = 0. Hence, C -1 will exist as noted earlier. The Det(C) is 
the square of a mesh area along a f = const plane. 

Cell Volume Specification 

The user has control of the grid by means of the initial surface point distri- 
bution and by specification of the cell volumes. Through the cell volumes 

the extent and clustering of the grid can be essentially controlled. Because the cell 
volume at each point must be given, it is clear that the user must devise some 
simple kind of method for specifying volumes There are many possibilities, and 
two approaches that have been used are illustrated here 

One way of generating mesh cell volumes is to form a grid about a ‘similiar’ 
but simple reference body for which the grid can be generated analytically, and to 
use the cell volumes from this reference grid for the more complex problem. For 
example, suppose we wish to grid an aircraft fuselage as a warped spherical-like 
grid following the topology shown in Fig.l To specify cell volumes we first find a 
sphere that has the same surface area as our fuselage and analytically build a grid. 
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For simplicity, choose a spherical grid that has uniform angle spacing and a radial 
distribution that varies exponentially In this special case, the control cell volumes 
are analytically known The grid cell volumes of this spherical reference grid are 
then used to specify the cell volumes of the fuselage grid. However, the fuselage will 
not have the same kind of surface area distribution as a sphere with equal angle 
distribution So here we use an adjustment of the form 


Al'i k, = 1(1 - v) - „| (Ai; , ,1,,*,,, 


(15) 


where u — > 1 for small l and u — + 0 for large / That is. the volumes are adjusted 
initially to the local boundary surface increments But as the solution is marched 
out, the u iform spherical volumes gradually become specified. Such an approach 
has been U'-ed, and. as a result, the far field portion of the grid tends to be uniformly 
spherical. Some of the results shown later will illustrate this behavior. 

Another method of specifying the volumes that has been used depends on 
the grid being generated and for this reason this second method is sometimes less 
robust. In this method, the specified volume at each point is set equal to the 
computed surface area element times a user specified arc length Specifically. 


where As is the user specified arc length and A 5 is the surface area. In all of 
our applications As has been given as an exponential stretching. In this kind of 
volume control specification, if an initial distribution of points is highly clustered 
in £ and 77, then these points tend to remain highly clustered in £ and 77 even far 
away from the body. To obtain a more uniform far-field distribution, the volumes 
specified from Eq. (16) have been averaged in £ and 77 with each step taken in g. 
For example 


AVj,Jt,i+i — (AVj+i^j-ri + AVj-ijfcjj.] - L AV’ :li ^ 1 j + i — A\' h k-i,i+i -f4AV’ J) fc ) j + i)/8 

(17) 

has been applied one or more times with each step in C- Alternatively, m the far- 
field the volumes given by Eq. (16) can be weighted with the uniform spherical 


volumes. 

Axis Treatment 

A coordinate singularity will be encountered whenever a whole grid face 
is mapped onto a closed body Here we will generate warped spherical grids, so 
Eqs. (3) become singular at the axis. With g the marching direction and using 
£ and 77 as sketched m Fig.l ( i.e. £ from pole to pole and 77 equatorial), then 
the axis at £ = 0 and the axis at £ = £ mai represent singularities. In particular, 
77-derivatives approach zero at the pole and Eqs. (3b) and (3c) are lost. There 
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are km ax points in the ^-direction. but. at a given $ -station, all x,y.z values are 
the same. Consequently, the difference equation corresponding to Eq. (3a) can 
be imposed kmax times to solve for onlv three unknowns, namely x, y. and z on 
the axis. Thus, the axis points are overdetermined and can be solved for via a 
generalized inverse scheme However, such a constraint is difficult to implement 
implicith into the solution algorithm 112] 

In the approach actually used at the axis, the bod\ is aligned so that y 
changes along the axis more rapidly than x or z We assume that the direction of 
the axis leaving the body surface f = 0 can be predicted from the surface distribution 
of points near the axis. For example, for the £ = £ m ax axis as sketched in Fig.2. 
averaging points in r/ at £ = £ max — Af and < or £ max — 2Af allows us to find an 
axis point inside the body from which the direction of the axis can be determined. 
Knowing the direction of the axis, changes of x and z along the axis can be expressed 
in terms of changes in y as 


x. = 


_ (V*) 
(Vy) 


V; 


Zr = 


11*1,. 

(Vy) 


(18a) 


(18b) 


where and are the known axis slopes and V implies backward differencing 
in the ^-direction. These relations along with the orthogonality relation provide 
enough information to back out the variables on the axis. Eliminating z- and z ? 
from the orthogonality relation Eq. (3a) by using Eq. (18) gives 


(Vx) 

.(Vy) 


6(X + 6^y + 


(If) 

(Vy) 


6{Z 


6<y = 0 


or since 6<y ^ 0 

+ = ° (19) 

Equations (18) and (19) can be readily built into the implicit solution process by 
modifying the end points of the f-block tridiagonal matrix formed by Eq. (12b). 
Let the axis points be denoted as j = 1 and j = jmax and define Mi — (Cf*Ai)/2. 
Then Eq. (12b) can be represented for a given index k from j = 1 to j = jmax as 
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( 20) 

where the contribution from the dissipation terms is not being shown The endpoint 
blocks F\. G’i. E }max and F jmax are supplied from the axis conditon, Eqs. (18) 
and (19), put into delta form as 


Vx 


{X],k,l+1 x J,k,l) l-l V],kl)—0 


k 1 4-1 - Zj k,l) - l - 1 - Vjkl)- -~^<{z 3 ,k,l - 1 - Zj,k.l ) = 0 


V A, 

— [V] k,l— 1 — y 3 ,k,l) ~ \Z) k l—l ~ z ] ,k,l) — 0 

v y 


Approximating Sc with Ac at j = 1. the blocks F\ and G i become 


F! = 


1 

Vi 

Vy 

0 


Vi 

Vy 

- 1 
_ Vi 

' v v 


0 

Vi 

Vy 

1 


(21a) 


Ci 


0 0 0 

V x ■» Vz 
Vv 1 V y 

0 0 0 


(21b) 


with g\ = 0. At a given f = const step, the inversion of a £-block tridiagonal 
results in slightly different predicted values of x, y, and z around the axis in rj 
Because these values should all be identical, their averaged value is used as the axis 
point. New axis slopes and are formed as each c = const grid surface is 
generated, and in this w r av the axis is allowed to curve. 

As a modification to what is described, second-order three-point differencing 
in £ is used (slightly weighted with the first-order differencing indicated above). 
Because three-point differencing is used at the ends, a modified block-tridiagonal 
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system that allows for additional blocks at locations (1.3) and (jmax,jmax - 2) is 
solved 

Boundary Surfaces 

If the grid is not being generated about a completely closed body surface as 
shown in Fig 1. then boundary conditons other than periodiciU or axis treatment 
must be imposed For example as sketched in Fig 3 it mat be necessary to 
generate a grid about a hall-span wing attached to fuselage or a plane of symmetry 
at. sat. £ = 0 In this situation, the grid at £ = 0 (see Fig 3i cannot be prescribed 
as this would be incompatible with the grid developing m the interior when marched 
in $ Consequently, the grid distribution on this boundary surface must evolve as 
part of the overall solution process There are multiple ways of generating the 
boundary grid along with the interior, one of which is to use one-sided differencing 
tied to the eigensystem of the C~ l A and C -i B coefficient matrices. 

The eigenvalues of C _ M and C~ } B have not been computed analytically, 
but. because of the special form of each coefficient matrix, one eigenvalue must be 
zero and the remaining two eigenvalues are real and of equal and opposite sign (see 
the Appendix). As a consequence, at least two proper combinations of the grid 
generation equations can be differenced to one side on a £ or rj boundary surface. 
Specification of the surface provides the third necessary condition to determine the 
grid points. At a £ = 0 boundary, for example, the characteristic combination 
of equations corresponding to the zero and the negative eigenvalues of C -1 A can 
be one- sided differenced. However, at this point we do not have analytical for- 
mula for the eigenvalues and eigenvectors of C -1 ,4 and C~ l B. and have, therefore, 
implemented another less rigorous approach. 

In lieu of working with the eigensystem. we have specified the surface and 
have simply discarded a governing equation. For example, the £ = 0 surface shown 
in Fig. 3 can be represented as 


f{x,y.z) = 0 (22) 

As long as this surface is close to a y = const plane. Eq. (22) can be used in place 
of the orthogonality condition. The remaining two conditions, restricted to a 
y = const plane, are essentially the two-dimensional governing equations, Eq. (2), 
so this approach would seem to a good one provided we can ensure compatibility 
with the interior mesh. To ensure this compatibility, we have built Eq. (22) into 
the solution algorithm, Eq. (11). as follows. Equation (22) is first differentiated 
with respect to f as 

f : ~ f x x ; f yVc f = 0 

and then approximating x- as x/ +1 — x/. etc gives 

f x (*e/ 4 i x i) fy(yi + 1 — yi) *r fz{zi- f*i — z i) = 0 (23) 
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Use of this relation in place of Eq (3a) leaves the coefficient matrix Bj and the 
vector gi of Eq. (9) unaltered, the first row of .4/ is replaced with zero elements, and 
the first row of C/ is replaced b\ the elements {f x .f y ,f z ) In the special case of a 
planar boundary, say y = 0. it can be verified that C7 1 .4/ has only zero eigenvalues. 
Therefore. 6c can be backward-differenced at £ = 0 without instability m this special 
case Backward differencing of dc at 5-0 for a more general surface f\x.y.z ) = 0 
with moderate curvature has also been stable 

RESULTS 

To demonstrate the capabilities of the current hyperbolic grid solver, grids 
were generated about two wdngs and about a simple wing-body configuration In 
all cases a spherical mesh topology was employed so an axis singularity exists. 

As a first illustration, a spherical-like grid was generated about a rectangular 
wing with a NACA 0012 airfoil section. The planform, the thickness distribution, 
and the airfoil section are shown in Fig. 4a where a nearly rectangular wing planform 
is defined from the superellipse 

(^/-t-mai) "f* (y/Vmax) = I 
and the thickness distribution is elliptical 

(y / Umax) ■+" (•Z/~maz)~ = 1 

Also shown is the user specified surface grid-point distribution which was arranged 
so that an axis extends from the back of the wing tips The surface grid has 79 points 
distributed in the ^-direction from pole to pole and 120 points in the periodic 77 - 
direction. The 5 = const surface lines were generated by computing the intersection 
of planes originating from a point located on the symmetry plane with the wing 
surface. The 77 = const surface lines were clustered as a cosine distribution. 

The grid was marched out in f using 40 steps to sweep out a distance of 
approximately 8 chords with a uniform grid spacing of 1/2 percent maximum chord 
specified as the first grid spacing off the body An entire warped spherical grid 
was generated without a plane of symmetry being used, but Figs.4b to 4e show 
only various segments of the generated grid for f and 77 = const planes. Portions 
of two different 77 = const planes in the z = 0 plane and in the vicinity of the 
axis are shown in Fig. 4b to indicate that the grid near the axis is relatively well 
behaved. The grid at the wing midspan ( y = 0) is shown in Fig. 4c. This is a 
5 = const plane. As seen, the grid is very smooth and uniform off the body and 
nearly orthogonal everywhere The tendency for the grids lines in the far-field of 
the grid to be circular reflects the fact that spherical mesh incremental volumes 
have been specified Finally. Figs. 4d and 4e show close-up views of this same grid 
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plane and illustrate that the grid generated near the profile nose and trailing edge 
regions is quite satisfactory 

A similar grid was computed as illustrated by Figs 5a to 5d in w’hich the 
grid spacing off the body is much finer and corresponds to the type of mesh spacing 
used in viscous flow simulation The w'ing m this case has a cambered airfoil section 
with the profile defined as 

- - 4~maz(\ x ~ x ) 

and a parabolic camber given b\ 

2 = 2 -*- 4cx( 1 - i) 

is sheared into the profile. Figures 5a to 5c show the plane of symmetry grid plane, 
w'hile Fig. 5d shows a £ = const plane in the vicinity of the the left wing tip. 

It must be remarked that the above wing calculation contained an airfoil 
with a sharp trailing edge, but unless the the airfoil trailing edge radius is rounded, 
the method can frequently break dowm. Provided the trailing edge radius is not 
zero, breakdown can be avoided by clustering a sufficient number of grid points in 
this region and adjusting the specified volumes to be compatible. 

Grids have also been generated for low aspect ratio wing-body combina- 
tions. and Fig. 6 shows portions of a mesh generated for a simplified space shuttle 
configuration which is missing the OMS pods, engine bells, and vertical fin. Figure 
6a show’s segments of the grid in the windward and leeward planes of symmetry 
(t] = const planes). A blowup of the nose region, showing the axis and canopy 
region, is provided by Fig. 6b. (The back of the shuttle was terminated in this 
case at a x = const plane.) Three cross sectional views (projected onto x = const 
planes) are provided by Figs. 6c-6d. To prevent grid line crossing in the corner of 
the upper wing body juncture region, it was necessary to run the grid solver with a 
value of 6 = 3 in Eq. (14). Flow calculations using these grids were recently carried 
out by Rizk and Ben-Shmuel [10]. 

CONCLUDING REMARKS 

A procedure has been developed for generating body fitted coordinates in 
three dimensions using hyperbolic partial differential equations. In this paper the 
scheme has been used to generate warped spherical-like grids about simple wing and 
wing-fuslage configurations. The hyperbolic grid generator can fail whenever the 
body surface is discontinuous or when the user specified surface grid distribution 
is too irregular. For simple continous body shapes, however, the hyperbolic grid 
generation scheme can be very fast and reliable. It requires a minimum of user 
interaction, and it can be used to generate grids suitable for either inviscid or 
viscous flow simulation. 
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APPENDIX 

Properties of the matrices C'~ l A and C~ ] B are needed to determine the 
hyperbolicity of the grid generation equations. These matrices can be written as 


{yr)$z ~ z ri l \y)y- ~ 1 (J/tjCz ~ 2 ri$y) 2 c + ?xs 

( 2 r>Cx ~ x n£z)y? ~ $y£yJ 1 ( 2 r)$x ~ x t t$z) z < "I" i£zJ 

( x r)$y ~ y^x)y- "t" Szsy*^ 1 [ x ri$y ~ Vr/Sx) 2 ' SztzJ 

and 

C~ l B 

( JDet C j) ” 

-r ^zC — yj?z)y> _r yysx*^ ' ( 2 *?y — y^sz) 2 ? yz$xJ 

( X S$Z - Ztfx) x ; T r)xSyJ~' (*({z ~ «{fa)y f + VySyJ ' 1 { x ‘Cz ~ z ‘$x) z s T r) z ^ v J ~ 1 

(y??z - x^y)x ( T [ytf x - Xc$y)y r - y y ? 2 J _1 (y<fx - *^y) 2 ? + riztzJ - 1 

where the Jacobian J = (AV) -1 is 

J~ l = xcy^z, x.ycZr, + x„y ? 2< - x^y f «„ - x^ycz, - x^y v zc 

and the transformation metrics which appear in the elements of A, B. and C are 

fx = J{yri z c — y$ z n) 

£y = J{ x ; z r) ~ x v z ‘) 

£z ~ J ( x nV ' — x cVt]) 
r\x = J{y- Z ‘ ~ Vi 2 ;) 
fly — J(xeZ- — 
r)z = J{x.yc - Xey.) 


C~ l A 

(JDet\C\) ~ 

(yr??z — z r)Sy) x c + fxsx*^ 
[ z rj$x ~ x r jCz) x ~ ~ Sy£xJ 
( X r/Cv yrjCx) x c 
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s* = J [yt z v ~ yv z d 

Si I — J{ x n~£ — * ‘~v) 

Cz = J(x*yr, - x„yc) 

Using the orthogonality relations and the definition the metrics sx-sy 
in terms of x.-.x,,. . it can be shown that the matrices C _1 .4 and C~ 1 B are 

symmetric For example, the (3.1) and (1.3) elements of C~ 1 A are found to be 
equivalent as shown b\ subtracting (3.1) from (1.3) lo obtain 

(VrjSz ~ Sis zJ ~ i ~ Vrjsx)^ ‘ xJ = 0? 

Using the definitions of J -1 £ x and J~ l £ z and eliminating terms which cancel leaves 

~~{ Z r) Z ' I" ““ - r r s’" ~ y?~rjfz = 0? 

Now using the metric definitions for g x etc 


[ZjrjZ^ X r) X-)(x 11 Ze XcZri) , (t/« Z^ Vr)Z‘ ) y$ s t} (^ Jl/ij Xrjye) — 0? 


Collecting the y^y* terms to complete the orthogonality relation • r- = 0 leaves 
two terms containing ye that sum to zero, that is 


~{zr,z. -t- y n y- - X, } x t ){x n ze - xeZrj) *r x v y,yez n - y^z v x n y^ = 0 


The remaining elements have also been checked for symmetry. The reader can also 
quickly verify that the degenerate case £ = x. y = y. and c — z gives the symmetric 
matrices 


C~ X A = 


0 0 
0 0 
1 0 


1 

0 

0 


and 


C~ X B = 


0 0 
0 0 
0 1 


0 

1 

0 


The matrices C~ l A and C _1 B have other interesting symmetries. In partic- 
ular. the determinant and the trace of of both matrices are zero. For example, the 
DetA is clearly zero from inspection of Eq. (6a) and Det(C~ l A) = DetADetC ~ l . 
The trace of C~ l A is given as 


trace(C A ) — {VtiSz z ri$y)X( -*■ frs xJ (Xrj^x "F 

Vri$x)Z( "F $z^zJ 
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or regrouping 


trace(C 1 A) = {y v ^z~z^ v )xr~{zr,<x- x n C 2 )y-^ (x n Cj / -y^Ci)2; + (?x^+f ! /$x-rb 2 ^ 2 ) J 

However. [z v y- - ynZ-)c 7 is -J~ } ixC z . etc . so the trace is zero 

Because 3 > 3 matrix C~~ 1 A is symmetric and has zero trace and determinant, 
one eigenvalue of C ~ 1 A must he zero and the other two real distinct eigenvalues 
must sum to zero Therefore, these nonzero eigemaiues must be of equal and 
opposite sign The matrix C~ l B has similar properties 
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a geometry and surface-point grid distribution on a rectangular wing 


Fig 4 Warped spherical grid about a rectangular wing 




b view of grid segments m the z =0 plane near the wing tip 


Continued 
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c. view of wing grid at midspan. 


Continued 
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d close-up of the leading edge midspan grid 


Continued 
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e close-up of the trailing edge midspan grid 


Continued 
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a. view of wing grid at midspan with cambered profile. 


Fig. 5 Hyperbolic grid generated about a rectangular wing with camber 
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b. close-up of the leading edge midspan grid 


Continued 
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c. close-up of the trailing edge midspan grid. 


Continued 
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1.00 



d. f = const plane near a wing tip axis. 


Continued 
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1000 



a. grid segments in the windward and leeward plane of the shuttle orbiter. 


Fig 6 Spherical-like grid generated for a simplified space shuttle orbiter. 
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b. Grid in the leeward plane near the nose of the shuttle orbiter. 


Continued 
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c Cross-sectional grid view of orbiter. near canopy, projected onto an x = const plane. 

Continued 
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d. Cross-sectional grid view of orbiter, strake region, projected onto an x = const plane 

Continued 
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e. Cross-sectional grid view of orbiter wing region, projected onto an x = const plane 


Concluded 
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